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ABSTRACT 

We present a new fast algorithm which allows the simulation of ionising radiation emit- 
ted from point sources to be included in high-resolution three-dimensional smoothed 
particle hydrodynamics simulations of star cluster formation. We employ a Stromgren 
volume technique in which we use the densities of particles near the line-of-sight be- 
tween the source and a given target particle to locate the ionisation front in the 
direction of the target. Along with one-dimensional tests, we present fully three- 
dimensional comparisons of our code with the three-dimensional Monte-Carlo radia- 
tive transfer code, MOCASSIN, and show that we achieve good agreement, even in the 
case of highly complex density fields. 

Key words: methods: numerical, radiative transfer 



1 INTRODUCTION 



Smoothed particle hydrodynamics (SPH) is a powerful and 
flexible numerical technique which has been applied to 
numerous astrophysical problems. The wide applicability 
of SPH derives largely from the fact that it is intrinsically 
Lagrangian. This allows SPH to deal effortlessly with large 
density contrasts and obviates the need for a grid or mesh 
, to be imposed on the simulation space, thus placing no 
restrictions on the symmetry of the problems which may 
be studied. SPH has b een used to m o del circumstellar 
and c ircumbinary disks llMurravj 119961). ISpeith fc Kunze 
20021'! ), stellar collisions (|Sills et all (|2002l h iLombardi et al 



,1996)), the collapse of m olecular cores (iBatd Il998l)). the 
forma tion of star clusters llBate et al. I (|2003T).|K1 essen et al.l 
l|2OO0h ) and galaxy formation ( jKobavashil l |2004 )). SPH is a 
proven tool for modelling gas dynamics, and hybridisation 
with N-body codes allows the inclusion of gravity in SPH 
simulations, but such comparatively basic input physics 
limit the sophistication of the simulations that can be 
performed using SPH. Two problems that have been 
traditionally very difficult to solve in SPH are the inclusion 
of magnetic fields and the implementation of radiative 
transfer. Studies of all the systems mentioned above are 
likely to benefit considerably from the ability to model 
these phenomena. 

In this paper, we describe a new algorithm designed to 
implement a very simple form of radiative transfer in the 
context of SPH simulations of the formation of star clusters. 
Massive stars inject energy and momentum into the ISM 



by a variety of means and these feedback mechanisms have 
important consequences on a variety of size scales, from 
the photoevaporation of circumstell ar disks around sing le 
or binary stars by ionising radiation l|Yorke fc Weld (|l996T )) 
to the inflation of galactic superbubbles by the combined 
action of the radiati on, winds and s upern ova explosions of 
thousand of O-stars (| Clarke fc Oevl (|2002h ). 

On a scale of a few parsecs, midway between the 
extremes given above, feedback from massive stars is 
thought to have a profound i mpact o n the evol u tion o f 
star clusters (e.g. Lada et all l| 19841 ) , lElmegreenl l|l998l ). 
iTenorio-Tagle et al.l ||2003|)) . Since most stars form in 
clusters I Clarke et al.l 1 2000i )). it is clear that the effect of 



stellar feedback acting on such scales is a crucial facet of 
the global star formation process. It has long been known 
that our galaxy is forming stars at only a fraction of the 
rate which would obtain i f star formation were a pro cess 
governed solely by gravity l|Zuckerman fc Evans! (Il974l )). It 
has also become clear that the majority of st ar clusters do 
not su rvive for long ti mes as shown b y (e.g.) lLada fc Ladal 
il2003l) in our Gala x y. iBastian et all (|2005h in M51 and 
iGoodwin fc Bastianl (|2006l ) in a sample of extragalactic 
young massive clusters. Feedback from massive stars offers 
attractive solutions to both of these problems. The ionising 
radiation and winds from massive stars drive powerful 
outflows into the gas around them. Such outflows can 
potentially remove enough mass from a protocluster on 
short timescales, leaving it unbound; the rapid expulsion of 
gas inhibits the production of stars, making star-formation 
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a very inefficient process. 

The influence of steiiar 
has b een stud ied from a 
files. i Goodwinl il 1997ft iHilL 



2001) 



feedback on star-formation 
va riety of theoretical an 



1980h iGever fc Burkeri 



Boil v fc Kroupa l|2003l ) and iGoodwin fc Bastian 



2006) investigated the effect of removal of the gaseous 



content of protoclusters, although they did not model the 
mechanism of gas expulsion itself, treating it either as 
instantaneous or as occurring on a prescribed timescale. 
Several studies have investigated t he dispersal of mol ecular 
clouds by photoionising radiation. IWhitworthl |l979') used 
simple theoretical arguments to study the erosion of a 
molec ular cloud by O-star s situa ted near the edge of the 
cloud, iTenorio-Tagle et all l|l986l ) used a one-dimensional 
Lagrangian finite-difference code to model the impact of 
photoionising radiation on the residu a l gas in globular 
clusters and iFranco fc Garcia-Segural |l996t ) conducted 
two-dimensional grid-based simulations of the growth of 
H II regions. All the theoretical work described above makes 
use of so me form of symmetry to si mplify the problem. More 
recently IGoodwin fc Whitworthl (|2004l ) adopted a fractal 
gas density distribution to mimic the asymmetries and 
inhomogeneities of star clusters. However, real protoclusters 
are inhomogeneous, anisotropic and dynamic environments. 
The only way in which to study ionising feedback under 
realistic conditions is to perform fully three-dimensional 
numerical simulations and SPH methods provide the ideal 
tool for this task. High resolution SPH simulations are 
computationally very expensive and the addition of further 
physics inevitably results in further overheads. However, 
with the advent of fast parallel m achines, such studies have 
now b ecome technically feasible. iKessel-Devnet fc Burkertl 
(2000) were the first to attempt this with their studies 
of radiation-driven implosion of molecular cloud cores 
|Kessel-Devnet fc Burkertl (|2003l )). 

In this paper, we present a fast and robust algorithm 
designed to simulate photoionisation in SPH and describe 
tests performed to ensure it produces accurate results. 
In Section 2, we give a brief overview of the physics of 
photoionisation. In Section 3 and 4, we describe how 
we have implemented the basic physics within the SPH 
formalism. In Section 5, we present one-dimensional tests 
of our algorithm and compare them to analytic solutions. In 
Section 6, we describe the comparison of the SPH algorithm 
to results generated with a fully three-dimensional Monte 
Carlo photoionisation code. Our conclusions are presented 
in Section 7. 



2 THE PHYSICS OF PHOTOIONISATION 

The idealised problem of the consequences of the sudden 
ignition of a source of ionising photons inside a cloud of 
uniform-density H I has been well studi ed, first by Stro mgren 
IStromgrenl jl939l 'l and later by Spitzer ISpiteeil |l978h . 

Initially, the ample supply of ionising photons causes an 
ionisation front (IF) to progress radially outward from the 
source at highly-supersonic speed leaving behind it a spher- 
ical H II region. During this phase of its evolution, the IF is 
said to be of R-type. Geometrical dilution of the radiation 
field and the recombinations of ions and electrons within the 



H II region progressively reduce the photon flux arriving at 
the front and slow its progress into the neutral gas. If the 
neutral atomic gas has number density no, ionisation of the 
gas will result in number densities of ions and electrons of 
rii and n e respectively. If the gas is taken to be H I clearly 
rii — n e and if the gas is fully ionised, then rii — n e = no. 
The recombination rate per unit volume is then given by 
an,n e = anj, where a is a recombination coefficient. If the 
source's ionising photon flux is Qh s _1 (where any photon 
whose energy exceeds 13.6 eV is regarded as ionising), the 
flux arriving at the ionisation front at radius Ri is given by 
Qh — F where F is the Stromgren integral given by 



F = 4tt 



r 2 n,QaBdr, 



(1) 



where no has been taken to be constant during the initial 
phase of the H II region's development and Ob, taken to 
be 3.0 x 10 -13 cm 3 s _1 , is the temperature-dependent 'case 
B' recombination coefficient which neglects recombinations 
directly to the hydrogen ground state. The ionising pho- 
tons produced by such recombinations are assumed to be 
absorbed elsewhere in the H II region and not to make any 
contribution to the ionisation equilibrium. This is commonly 
referred to as the 'on-the-spot' (OTS) approximation and is 
justifiable in cases where the optical depth of the H II region 
to secondary ionising photons is smaller than the dim ensions 
of the H II region. Yorke in iKudritzki et al. I (|l988l ) shows 
that, for an O-star with an effective temperature of 40000 
K, the on-the-spot approximation is valid throughout most 
of the H II region at the densities of interest here (> 10 2 
cm" 3 ). 

As the ionisation front proceeds into the neutral gas 
surrounding the source, the point is quickly reached where 
Qh = F and the ionisation front can proceed no further. 
The radius at which this happens is known as the Stromgren 
radius, Rs and is defined by integration of Equation 1: 



Rs = 



3Q H 

\--KTi\olb ' 



(2) 



Since the ionised gas within the H II region is expected 
to reach an equilibrium temperature of ~ 10 4 K, whereas 
the ambient temperature of the neutral hydrogen in a typ- 
ical GMC is ~ 10 K, a very large pressure gradient then 
exists across the ionisation front. This will cause the H II 
region to expand violently into the surrounding H I, driving 
a strong shock before it. This expansion is negligible during 
most of the initial phase described above, since the speed 
with which the ionisation front propagates is much greater 
than the speed of sound in the ionised gas (~ 10 km s _1 ). 
However, once the ionisation front has slowed down to a 
speed of a few times the sound speed, the pressure- driven 
expansion of the H II region begins to dominate the evolu- 
tion of the system. The IF is then said to be D-type. 

By considering the jump cond itions across b oth the 
shock front and the ionisation front, ISpitzerl (|1978l ) derived 
a simple expression for the time-dependence of the radius of 
the ionisation front, Ri, during the expansion phase of the 
H II region 



Ri{t) = 



Rs ( 



1 + l£fl\ 
ARsJ 



(3) 
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where c 3 is the isothermal sound speed within the H n 
region. 



3 BASIC STROMGREN VOLUME METHOD 

We have conducted our simulations using a sm oothed par- 
ticle hydrodynamics code, fully described in [Bate et al.l 
l| 19951 ). The code is a hybrid N-body/hydrodynamic code, 
using a binary tree algorithm for gravitational calculations, 
and the SPH formulation to handle the gas dynamics. Par- 
ticles are evolved on individual timesteps to improve effi- 
ciency. This code is a proven tool fo r study i ng clustered star 



forma ti on (see, e.g 



proven tool tor studying clustered star 
Bon nell et ail (|200lh , iBonnell fc Bate! 



(|2002h . [Bate i et all (|2003h ) and is thus ideal for studying the 
problems outlined in the Introduction. 

Stars are represented in the code by sink particles. A 
sink particle is a point mass with an accretion radius. Gas 
particles straying inside a sink particle's accretion radius are 
accreted by the sink particle if they pass a series of the tests, 
the most obvious of which is whether or not the gas parti- 
cle is gravitationally bound to the sink particle. If the gas 
particle passes the test, its mass, momenta and energy are 
added to those of the sink particle and the gas particle is 
removed from the simulation. 

In previous simulations, the sink particles interact with 
the gas particles only via their gravitational fields and the 
accretion process described above. We have modified the 
code to treat one or more of the sink particles within a 
simulation as a source of radiation. The photoionisation al- 
gorithm selects a gas particle as a target and determines 
whether or not enough ionising photons from a given ra- 
diation source reach the target particle during the current 
time-step to ionise it (or, if it is already ionised, to keep 
it that way). Consistent with the physical assumptions de- 
scribed above, all gas particles in the simulation are regarded 
as being fully ionised or fully neutral. Once the ionisation 
algorithm has decided which particles are ionised, these par- 
ticles are heated to a temperature of ~ 10 4 K appropriate 
for the choice of as, and control is returned to the dynam- 
ical portion of the code. 

In a nebula consisting of pure hydrogen and charac- 
terised by conditions where the on-the-spot approximation 
is valid, secondary ionising photons do not contribute to the 
overall ionisation equilibrium. The photons which do con- 
trol the ionisation balance are therefore the primary ionising 
photons which travel radially outwards from the radiation 
source until they are absorbed. Since these photons follow 
purely radial trajectories through the gas, the fate of each 
photon is independent of those traveling along other radial 
trajectories. This has the consequence that, in a non-uniform 
cloud, each photon can be treated as though it were moving 
through a spherically-symmetric cloud possessing whatever 
radial density profile exists along the photon's direction of 
propagation. If the cloud in question is everywhere suffi- 
ciently dense that the on-the-spot approximation holds, the 
ionisation structure around a point source of radiation can 
then be found by a Stromgren volume method in which the 
distance from the source to the ionisation front in any given 
direction can be found using the radial density profile in 
that direction. It is worth noting at this point that in the 
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Figure 1. Illustration of the method used to select particles 
whose densities are to be used to derive the density profile along 
the line-of-sight between the radiation source and a given target 
particle, denoted by the thick dashed line. 



case of multiple ionising sources with overlapping Stromgren 
spheres, the algorithm described above cannot reproduce the 
radiation field in the overlap region. The ionised mass will 
therefore be underestimated by the current methods. This 
may have important consequences in the case of crowded 
fields and therefore a new version of the photoionisation al- 
gorithm is currently being developed and will be described 
in a forthcoming paper. 

In a spherical cloud with an arbitrary radial density 
profile p(r), Equation 1 becomes 



F = 



4irr 2 n(r) 2 ceBdr 



(4) 



where the uniform number density no has been replaced with 
a number density n(r) which is a function of radius. 

Hence, the flux arriving at any give particle can be 
found if the function n (r) along the line-of-sight t o tha t par- 
ticle can be estimated. iKessel-Devnet fc Burkertl (|2000t ) pro- 
posed an ingenious method by which the function n(r) can 
be evaluated in SPH and we make use of it here, with some 
modifications. The method is illustrated in Figure [T] A tar- 
get SPH particle is first selected and a line-of-sight from this 
target to the source is drawn. The target's neighbours (those 
inside the dotted circle in Figure [1] are fetched from the SPH 
code's neighbour-lists and the one closest to the line-of-sight 
is selected (labeled as number 4 in our diagram), that is, the 
neighbour for which the angle d marked in Figure [J is the 
smalle st. To save computer time. IKessel-Devnet fc Burkertl 
(2000) discussed at length the use of a tolerance angle O, 
such that the first particle in a given target's neighbour list 
which satisfies 9 < O is selected, but we found that the per- 
formance of our algorithm remains acceptable even when it 
is forced to search for the closest particle to the line of sight. 
The position of the selected particle is then projected onto 
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the lin e-of-sight to produce what iKessel-Devnet fc Burkertl 
(2000) term an evaluation point, labelled r<t in Figure [1] The 
neighbours of the newly-chosen particle are then fetched 
and the selection process is repeated until the radiation 
source is reached. Once t his pr ocess has been completed, 
IKessel-Devnet fc Burkertl (|2000h then solve an equation for 
the rate of ionisation of the target particle. Instead, we use 
the density profile determined by the process outlined above 
in a discretised version of Equation!?] to decide whether the 
target particle receives enough radiation flux to ionise it, 
a faster and more robust meth od than that employed by 
IKessel-Devnet fc Burkertl (|2000t ). 

From the densities of the Nlos particles selected along 
the line-of-sight our target, we know the value of n(r) in a 
series of Nlos radial bins whose inner and outer radii are 
defined by the projected positions of the particles along the 
line. We assume that the densities at these projected posi- 
tions are simply the same as the densities of the relevant 
particle i, evaluated in the usual way as 

N ni!igh 

Pi = m,iW(0,hi)+ 2J mjW(rij,hij), (5) 
i=i 

where m denotes the mass of each particle and the sum runs 
over the N neig h neighbours of particle i (defined as all other 
particles within two smoothing lengths, 2hi, of particle i). 
W is the smoothing kernel, a function of the interparticle 
separations rij and the average smoothing lengths Tiy of 
particle i and each of its neighbours. 

Alternatively, we could estimate the densities at the 
evaluation points by performing an SPH density estimation 
of the form 

JV 

Pk = y^m j W(r jk ,h j ), (6) 
i=l 

where the sum instead runs over all N SPH particles that 
overlap the fe-th evaluation point, rjk is the distance of each 
of these particles from the evaluation point and hj is the in- 
dividual smoothing length of each particle. However, iden- 
tifying the N particles is not trivial; some but not all are 
likely to be neighbours of the fc-th particle on the line-of- 
sight to which the evaluation point belongs. Finding them 
is therefore time-consuming and we find that in any case 
it makes very little difference to the results; the density of 
the fc-th particle is almost always a good estimator of the 
density at the fc— th evaluation point. 

We now approximate the integral in Equation [4] by a 
sum as follows: 

— = rUn(n)) 2 a B An. (7) 

i=l 

The inner radius of bin i is defined by the radial position 
of particle i — 1 projected onto the line-of-sight, r(i — 1), 
and the outer radius of bin i by the projected radial 
position of particle i, r(i). The width of each bin is then 
r(i) — r(i — 1) = Ar^. (n(n)) is the average number density 
of the two SPH particles whose positions define the bin, 
[n(i) + n(i — l)]/2. The density of the point mass, used 
in calculating the number density of the first bin, is taken 
to be zero. Once the sum has been calculated, F can be 
subtracted from the source luminosity Qh to obtain the 



flux of photons arriving at the target particle. We can then 
determine the particle's ionisation state during the current 
time- step. 



4 IMPROVEMENTS OVER THE BASIC 
STROMGREN VOLUME METHOD 

Although the method described above is adequate for the 
simulation of the growth of a spherical H II region in a 
uniform motionless cloud, our code is also intended for 
use in highly-inhomogeneous and dynamic situations. It 
is highly likely under such circumstances that neutral gas 
may enter the H II region from outside, in an accretion flow, 
for example, and that ionised gas may either leave the H II 
region, or be otherwise cut off from the supply of photons 
required to prevent it recombining. We have therefore 
improved upon the basic Stromgren volume method to cope 
with these eventualities. 



4.1 Neutral gas inside the H n region 

Modifying the method above to account for the photons re- 
quired to ionise neutral material which has entered the H II 
region (perhaps via an accretion flow), or to grow the H II 
region from scratch (the R-type phase of the IF's evolu- 
tion) is not difficult. If we introduce an ionisation fraction 
Xi (still taken to be either or I) for each particle, then we 
can rewrite Equation [7] as 

N LO S N L O S . 

±= rnn(n)) 2 a B Ar l+ ]T r?(n(r i ))(l - a*)^. (8) 

i=l i=l 

The first term in the equation, the flux required to keep pace 
with the recombinations in ionised gas, is always present. 
The second term allows for the possibility that some flux 
may be consumed by ionising the gas from a neutral state 
first, and makes the expression time-dependent. The term 
An /At is the discretized speed at which the ionisation 
front propagates. The timestep on which the algorithm 
evolves the ionisation front, At, is usually set to the shortest 
dynamical timestep in use by the SPH code at the time, so 
that the H n region is updated whenever any SPH particles 
are updated. If this timestep is very short, the flux arriving 
at a given target particle may not be sufficient to ionise 
it immediately, in which case the flux is banked - allowed 
to accumulate until the particle can be ionised. Particles 
with banked flux are not heated above the neutral gas 
temperature. 



4.2 Ionised gas straying outside the H II region 

Ionised material that is for some reason deprived of photons 
will recombine and eventually cool. To allow for this 
eventuality, ionised particles which are receiving no ionising 
flux are flagged and monitored by a neutralisation routine. 
This routine reduces the ionisation fraction of such particles 
according to their local recombination timescale, (an) -1 . If 
the ionisation fraction of a particle drops below a half, it is 
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Figure 2. Decision tree for the SPH ionisation code. 



ass umed to become fully neut ral. 

iKessel-Devnet & Burkeru (|200Gh note that cooling 
recently-recombined gas instantaneously back down to the 
original ambient gas temperature is not a good treatment 
of the physics of recombination, since recombined gas, al- 
though neutral, is likely to still be hot. Neutralised particles 
are therefore passed on to a cooling routine which assigns 
them new t emperatures according to t he co oling curve 
described by ISchmutzler fc Tscharnuterl (ll993T l until they 
reach a temperature of 100 K, at which point they are re- 
turned to the initial average gas temperature (usually 10 K) . 



4.3 The ionisation code's decision tree 

Once the flux reaching a given target particle has been deter- 
mined, taking into account the consumption of photons by 
neutral material lying on the line-of-sight from the source, 
the ionisation state of the target particle is determined us- 
ing the decision tree illustrated in Figure[2] Note that in this 
tree, the approximation is made that if the flux reaching an 
ionised particle is non-zero, that particle remains ionised. 
It could be argued that the flux must in fact be sufficient to 
exceed the recombination rate in the particle. However, in 
all cases, the fluxes calculated by the code are those reaching 
the centres of particles and the decrement in flux incurred 
between the particle edge and its centre is taken account of. 
The only particles for which this assumption is questionable 
are those at the ionisation front. For a given ionised particle 
i at the ionisation front, a non-zero flux reaches the centre 
of particle i but zero or 'negative' flux reaches the next par- 
ticle further out along the line-of-sight to the source, j. This 
means that the true ionisation front lies somewhere between 
the centres of i and j and thus somewhere inside both par- 
ticles (recall that SPH particles overlap one another). The 
decision to keep i ionised and j neutral introduces an error 
in the location of the ionisation front along the line-of-sight 
on which i and j lie of ~ hi , consistent with the accuracy to 
which any spatial quantity can be meaningfully determined 
in SPH. 
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Figure 3. Evolution of the radius of the ionisation front with 
no heating of ionised particles, comparing the approach to the 
Stromgren radius produced by the SPH ionisation code (dashed 
line) to the analytical solution given in Equation 1101 (solid line). 
The Stromgren radius itself is shown as a dotted line. 



5 ONE DIMENSIONAL TESTS 

The most basic tests one can perform with our hybrid hy- 
drodynamics/ionisation code are to follow the approach to 
the Stromgren radius and the subsequent expansion of the 
H II region. 

The sudden ignition of a source of Qh ionising photons 
per second in a uniform cloud of atomic number density 
no leads to the highly supersonic expansion of an R-type 
ionisation front. The expansion velocity of the front can be 
derived by considering the rate of fresh ionisations at the 
front. If the IF has an instantaneous radius Ri and the flux 
reaching the front is Ft, we may write 

dJ?7 Fj 1 



4 2 n 2 

-7rn itjaB 



(9) 



This can easily be integrated and combined with Equation 
rjto yield 



Ri(t) =R S (1- exp [-n a B t]) 3 



(10) 



In Figure we show the results of the ignition of an 
ionising source in a uniform cloud, following the approach 
of the ionisation front to the Stromgren radius. To simulate 
this, the heating of ionised particles was disabled and the 
timestep on which the code dumps output was set to a very 
small value. We see that the agreement between the SPH 
results and the analytic solution is excellent. 

As the ionisation front approaches the Stromgren 
radius, its velocity slows towards the speed of sound in the 
neutral material surrounding the H II region. The IF then 
transitions to its D-type phase according to Equation [3] 
|Spitzerl |l978)). To test the ability of our code to replicate 
the Spitzer solution at different numerical resolutions, we 
constructed three spherical clouds of identical mass M c id 
and radius R c id, containing ionising sources of identical 
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Figure 4. Evolution of the radius of the ionisation front in the 
LR (dotted line), MR (dashed line) and HR (dash— dot line) cal- 
culations, compared with the Spitzer solution (solid line). 



luminosities Qn, but built form different numbers of 
particles. Since the Stromgren radius in the three clouds 
is the same, the initial H II regions in the clouds contain 
different numbers of particles. We chose M c id = 17OOM0, 
Raid = Ipc, Qh = 10 49 s-\ so that n = 1.65 x 10 4 cm" 3 
and Rs = O.lpc. We chose particle numbers of 10 6 (run 
HR), 10 5 (run MR) and 10 4 (run LR), so that the initial 
Stromgren radii contain respectively 1000, 100 and 10 
particles. 

It is not a priori obvious what the minimum re- 
quirement to resolve an H n region in SPH might be. To 
resolve an object in SPH is usually taken to require at 
least the average number of particles c ontained within 
the S PH smoothing kernel (usually 50). iBate fc Burkertl 
1 19971 ) found that, to resolve gravitational fragmenta- 
tion correctly required that a Jeans mass contain at 
least 100 particles (twic e the canonical resolution limit). 
However, iHubber et all l|200rj) find that even seriously 
under-resolving the Jeans mass (with < 10 particles) 
does not result in spurious fragmentation. Our run LR 
also under-resolves the H II region, in the sense that ten 
particles is well below the canonical SPH resolution limit, 
meaning that the H II region may behave unpredictably. 

In Figure [4] we plot the time-evolution of the ioni- 
sation front radius in our three models together with the 
Spitzer solution for comparison. We see that the three 
SPH models agree very well both with each other and 
with the Spitzer solution, at least for the evolution of the 
H n region from Rs to ~ 10i?s- Even the low-resolution 
calculation gives reasonable agreement with the analytical 
solution over the range of radii considered. For this simple 
problem at least, our implementation can faithfully follow 
the evolution of an H II region which would be considered 



severely under-resolved by the usual SPH resolution criteria. 



6 COMPARISON WITH A 3D MONTE CARLO 
PHOTOIONISATION CODE 

Since we intend to use our ionisation code in situations 
where the gas density distribution is highly inhomogeneous, 
one-dimensional tests can only take us so far. In order 
to perform a more stringent test than those presented in 
the previous section, we therefore selected a very complex 
density distribution derived from a three-dimensional 
simulation of clustered star formation, and compared the 
H II regions generated by our code and by the fully 3D 
Monte Car lo photoionisation and rad iative transfer code, 
mocassin (|Ercolano et al.1 [20031 . 120051 ). This code adopts a 
stochastic approach to the transfer of radiation, by simulat- 
ing the individual processes of absorption and re-emission 
of radiation quanta as they leave an ionisation source and 
diffuse through a gas. The polychromatic nature of both 
the primary (stellar) and secondary (diffuse) components of 
the radiation field are accounted for self-consistently as the 
photon trajectories are computed and the ionisation and 
thermal structure of the nebula established. All relevant 
atomic physics processes are included and the atomic 
database is updated reg ularly, currently u sing a combi- 
nation the Chianti V5 (|Landi et al.l 2006) database for 
collision strengths, transition probabilities and energy levels 
of heavy metals, the opacity project for photoionisation 
cross-section and a new calculation of the H + , He + and 
He 2+ free-bound continua presented bv lErcolano fc Storevl 

Dust grains mixed with gas inside the ionised regions 
may also affect the global ionisation structure of the gas. 
The grains compete with the gas for the absorption of UV 
radiation, hence reducing the flux available for the pho- 
toionisation of atom and ions. This can be self - consis tently 
accounted for by mocassin (1 Ercolano et al. | l2005h . but 
not by the SPH ionisation algorithm, since the amount of 
radiation subtracted from the field, and thus available to 
the gas, depends heavily on a number of variables, including 
the local dust-to-gas ratio, the grain chemistry and size 
distribution, as well as the shape of the local radiation 
field, which cannot be obtained by the SPH ionisation 
algorithm. In order to be able to compare the outputs from 
the two codes we did not include any dust in the MOCASSIN 
calculation, nevertheless we point out that the ionised 
masses thus obtained may be slightly overestimated. 

We chose to compare MOCASSIN and the SPH ionisation 
algorithm under realistic astrophysical conditions. The gas 
distribution selected for the comparison was drawn from an 
SPH simulation of t he collapse of a molecu lar cloud orig- 
in ally performed by iBonnell fc Bate! l|2002l ) and repeated 
bv lDale et ail (2005), who modeled the effect on the cloud 
and the stellar cluster formed by it of the ionising radiation 
from the single O-star born at the cluster centre (using the 
code described here). Three views of the gas distribution 
are shown in the top panels of Figure [6] The ionising 
source is at the centre of the image, at the intersection 
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of several dense filaments of gas. The density structure of 
the cloud is evidently highly inhomogeneous - maximum 
and minimum densities are ~ 2 x 10 7 and ~ 2 x 10 1 cm" 3 
respectively, with mean and median densities of ~ 2 x 10 4 
and ~ 4 x 10 2 cm -3 respectively (these densities are derived 
from the densities of the SPH particles). Note that, as well 
as exhibiting a very large dynamic range, these densities 
are sufficiently high that the OTS approximation will hold. 
Modelling photoionisation in this highly inhomogeneous 
(and, from the source's point of view, highly anisotropic) 
environment is clearly very challenging. 

MOCASSIN is a grid-based code, in which a 'mother- 
grid' of arbitrary resolution is defined; any cell in this grid 
may be further divided into a finer subgrid of arbitrary 
resolution, so that arbitrary gas distributions can be effi- 
ciently represented in the code with adequate detail. SPH 
is a gridless method, representing gas distributions by an 
ensemble of overlapping particles of (usually) equal mass, 
but varying radius (the smoothing length). We therefore 
devised a method to interpolate the SPH gas distribution 
onto a grid suitable for use by MOCASSIN and subdivided in 
such a way that the resolution of the grid was everywhere 
as close as possible to the resolution of the SPH particle 
ensemble. We first generated a uniformly spaced 63 3 mother 
grid with cell size d max filling the same volume as the SPH 
particle distribution. Each SPH particle was then located 
in the grid and the smoothing length of the smallest 
particle, h m i n , (if any) in each mother grid cell determined. 
If for any cell d max > 2h m i n , the cell is subdivided into 
{1 + INT[d maa; /(2fo m i n )]} 3 cells. Finer subdivision of grid 
cells, for example by h m in or h m i n /2 did not result in a 
significant improvement in the interpolation accuracy. This 
procedure resulted in the division of 223 of the 250 047 
mother grid cells into subgrids. Once the simulation volume 
is suitably partitioned, the mass of each cell is calculated by 
finding all N ce a SPH particles overlapping the centre of the 
cell, estimating the density at that point using Equation [6] 
with the sum running from 1 to N ce u, and multiplying the 
density estimate by the cell volume. 

A simple check reveals that the total mass in our SPH 
particle distribution is 514.8 Mq, compared to 513.3 Mq in 
the adaptive grid, implying that our interpolation method is 
adequate. In Figure [6] we compare column densities maps, 
as seen along the three principle axes, of the SPH particle 
distributions (top row) and the density grids created for 
MOCASSIN (bottom row). We see that the MOCASSIN grids 
appear to reproduce the morphology of the SPH particle 
distribution very well. It should be noted that, in order to 
facilitate visualisation, the information contained by each 
subgrid and SPH particle has been averaged and mapped 
onto a uniform pixel map, so that the resolution of the 
figures shown is much lower than that actually used in our 
calculations. 

To compare the H II regions generated by the two 
codes, we placed sources emitting Qn = 10 49 (ionising) 
photons s _1 at position (0, 0, 0) in both simulation volumes. 
The modifications described in Section 4 were disbaled, 
so that the SPH code was operating in a pure Stromgren 
volume mode. We first used extremely metal-poor abun- 
dances for the MOCASSIN run. The general H II region 
morphologies predicted by the two codes are shown in 
Figure [7] and are in good agreement, although a closer look 



will reveal low density ionised protrusions in the SPH maps 
that are absent or much less accentuated in the MOCASSIN 
maps. In terms of total ionised mass fractions predicted, 
however, the morphological discrepancies seem to have 
little significance, with total ionised masses being predicted 
by the SPH algorithm and the MOCASSIN code of 20.9M Q 
and 2O.5M0 respectively. The extra features appearing in 
the SPH maps in fact tend to be of much lower densities. 
The qualitative morphological differences but quantita- 
tive agreement may be caused by the inherently poorer 
resolution of low-density regions by the SPH code, or it 
may be an artifact of the interpolation of the Lagrangian 
SPH mass distribution onto a set of Eulerian grids. How- 
ever it is more probable instead that temperature effects 
may be the dominant cause, as described in the next section. 



6.1 Temperature-dependance of the 
recombination coefficient 

The isothermal nature of the SPH ionisation algorithm in- 
troduces a certain level of uncertainty in the calculation 
of the ionised region due to the temperature dependence 
of the recombination coefficient. The MOCASSIN model pre- 
sented here, where the gas temperatures are calculated self- 
consistently with the ionisation structure, shows that the 
ionised region of such a complex density field is far from 
being isothermal. The temperature structure of the ionised 
region is illustrated in Figure [5] as the histogram of mass 
fraction (black solid line) and the cumulative mass fraction 
(red dashed line) of material at a given temperature. It is 
clear that the success of the SPH algorithm to calculate the 
correct ionised mass of the grid relies on the careful choice 
of the value of the recombination coefficient, and implied 
gas temperature. For a density inhomogeneous pure hydro- 
gen nebula of arbitrary geometry, the total ionised volume 
is obtained from the solution of 



N p N e a(T)-dV = 



(11) 



where, for almost complete ionisation, N p —N c ~Nn, the 
hydrogen number density. The total ionised mass of a region 
under such assumptions, therefore, strongly depends on the 
recombination rate averaged temperature, defined as, 

f v n h T-dV 



< T[Njj] >-- 



I v N 2 H -dV 



(12) 



The MOCASSIN calculations yield < T[N%] > = 9120 K, 
which is very close to the value 8950 K implied by the choice 
of olb = 3.0-10~ 13 cm 3 s _1 in the SPH calculation. The 
proximity of these two temperatures allowed for the excellent 
agreement in the global ionised masses calculated by the 
two codes. As anticipated in the previous section, however, 
the detailed ionised volume shapes do show some differences 
suggesting that the temperature structure along a given line 
of sight may deviate from that implied in the SPH. 

The temperature at any point in the ionised region is 
determined by the balance between the heating, dominated 
by the mean energy of the photoelectrons absorbed by H 
and He, and cooling, which when metals are present, is 
dominated by collisionally excited emission lines (particu- 
larly infra-red fine-structure lines of [O ill]). In the absence 
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Gas temperature T [kK] 



Figure 5. Mass fraction (black solid line) and cumulative mass 
fraction (red dashed line) of material at temperature T in the 
ionised region. 



of metals, the cooling is dominated by collisionally excited 
Lya, which is much less efficient at removing energy from 
the system, therefore metal-rich regions show lower overall 
temperatures and steeper temperature gradients than their 
pure-H counterparts. Metallicity will therefore influence 
the temperature structure, and consequently the ionised 
mass of a region. This can be taken into account in the 
SPH ionisation algorithm by choosing a different alpha 
value appropriate for each case. A good estimate of the 
alpha value can be obtained by running snapshots of the 
evolving density structure calculated by the SPH code at 
some representative ages through a 3D photoionisation 
code, such as MOCASSIN, to self-consistently calculate the 
temperature and ionisation structure and thus obtain the 
recombination rate average temperature for the grid. The 
appropriate alpha value for the SPH ionisation code will 
therefore be that calculated at the recombination rate 
averaged temperature found. 

We should finally note that the range of gas temperature 
obtained throughout the ionised region may indeed affect 
the hydrodynamical evolution of the system (the SPH 
ionisation algorithm assigns the same gas temperature to 
all ionised particles). We defer the investigation of these 
potentially important effects to future work. 



6.2 Effects of dust 

None of the calculations presented here include the effects 
of dust, which would in reality also act as a s ink of ionising 
photons - for example, IWood fc Churchwelll l|l989h suggest 
that up to 90% of the stellar flux in an HII region may be ab- 
sorbed by dust. It is not clear what effect dust would have on 
the morphology of complex HII regions such as those studied 
here. If the dust distribution closely follows that of the gas, 
the effect would essentially be to shrink the HII region but 
to preserve its shape. This effect could be be approximately 
modelled by modifiying the recombination coefficients used 
in the codes. However, if the dust and gas are not strongly 
correlated, two different radiative transfer problems must be 
solved simultaneously. In this the coupling of the radiative 
transfer solutions to the hydrodynamic evolution of the sys- 



tem in question is much more complicated, but we do not 
attempt to address this problem here. 



7 CONCLUSIONS 

We presented the results of tests of a new fast algorithm 
for simulating ionising radiation from point sources in the 
context of smoothed particle hydrodynamics (SPH) simula- 
tions. The method we use is essentially a Stromgren volume 
technique. We use a met h od si milar to that presented in 
iKessel-Devnet fc BurkerU (|2000l ) to estimate the density 
profiles along lines of sight to individual SPH particles 
before performing a Stromgren integral to calculate how 
much flux, if any, reaches that particle. We improved on the 
simple Stromgren volume method by including in the flux 
calculation the flux required to ionise any neutral material 
lying on the path to the target particle. This modification 
is relevant in problems were the dynamical timescales can 
be as short as the timescale required to ionise the gas, 
which can be the case in accretion flows. We also developed 
a method to allow ionised particles which are somehow 
deprived of their photon supply to recombine and cool using 
an optically-thin cooling curve. This modification is also 
important in dynamical applications where gas can move 
out of an H II region, or be shadowed by dense material 
intervening between it and the radiation source. 

We conducted simple one-dimensional tests in which 
we showed that the algorithm was able to reproduce the 
ionisation front's R-type phase, in which it approaches the 
Stromgren radius at very high velocity, and the D-type 
phase in which the H II region expands thermally at 
supersonic but ever-decreasing speed. We showed that the 
SPH code is able to reproduce the D-type expansion phase 
well even when the number of particles inside the H II 
region is well below the canonical SPH resolution limit. 

Having demonstrated that the algorithm is able to 
reproduce the analytical solutions to simple problems, 
we compared the results generated in a highly inhomo- 
geneous and anisotropic gas distribution to the output 
from MOCASSIN, a much more sophisticated Monte Carlo 
photoionisation and radiative transfer code. We found that, 
for an appropriate choice of the recombination coefficient 
a, the two codes agreed very well, to within 2 percent, 
on the quantity of gas that should be ionised, and agreed 
reasonably well on the morphology of the H II region. 
We conclude that the algorithm is adequate for modelling 
photoionisation in SPH simulations of star formation with 
very complex gas distributions. 
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